Optimizing Replica Exchange Moves For Molecular Dynamics 



Walter Nadler^E and Ulrich H. E. Hansmann 2 - 

1 John-von-Neumann Institute for Computing, Forschungszentrum Julich, D-52425 Julich, Germany 
2 Department of Physics, Michigan Technological University, Houghton, Michigan, USA 

(Dated: February 1, 2008) 

In this short note we sketch the statistical physics framework of the replica exchange technique 
when applied to molecular dynamics simulations. In particular, we draw attention to generalized 
move sets that allow a variety of optimizations as well as new applications of the method. 
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Effective simulation of proteins [lj, glasses @, and 
similar complex systems [3] is hampered by slow relax- 
ation due to barriers and bottlenecks. Replica exchange 
0, IE Q - also known as parallel tempering - is one of the 
main approaches to overcome these problems. Originally 
devised for stochastic simulations, it is used nowadays 
also in combination with molecular dynamics (MD) simu- 
lations, i. e. simulations that have a strong deterministic 
character. We have found that there exists some confu- 
sion on the correct application of the replica exchange 
approach to MD. In this short note we will sketch the 
theoretical basis for applying replica exchange to MD 
simulations, introduce generalized move sets necessary 
for optimizing exchange rates, and point to possible ex- 
tensions of the concept. 

In replica exchange, a set of stochastic simulations 
are performed in parallel with distinct weight functions 
Wi(S), see Fig. 1. At certain times an exchange of current 
conformations Si of replicas at neighboring samplers Q 
is attempted, e.g. for a single pair of simulations 
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FIG. 1: Sketch of the basic situation in replica exchange sim- 
ulations: simulations run independently on stochastic sam- 
plers with weight function Wi(S); at certain times exchange 
of the current conformations is attempted with probability 
p M , Eq. ©. 



{Si, S2} — > {Si, S' 2 } = {52, Si} 



(1) 



Such an exchange move does not change the statistics of 
the full distribution function 



w to t{Si,S 2 ) = wi(Si)w 2 (S 2 ) 



(2) 



if it is accepted or rejected according to a generalized 
Metropolis rule [§], 



PM ({5i,5 2 } -» {S[,S 2 }) = min 



w t Qt{S'i, S' 2 ) 
' wtot(Si,S 2 ) 



(3) 



In physical and chemical applications one usually fo- 
cuses on the canonical ensemble. Using replica exchange 
an individual replica performs a random walk in temper- 
ature, allowing it to enter and escape local potential min- 
ima. As a consequence, the state space is explored more 
thoroughly, especially at low temperatures. The weight 
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functions employed in such situations are the canonical 
distributions, w(S) oc exp [— f3E(S)], with E(S) the en- 
ergy of the system and [3 = X/ksT the inverse temper- 
ature. For the exchange move Eq. (JT|), the Metropolis 
criterion assumes the simple form 

PM {{Si, S 2 } {5 2 , 5i}) = min (1, exp [A/3AE]) (4) 

with A/3 = fa - fa and AE = E{S 2 ) - E{Si). 

The exchange move of Eq. |TJ is not unique. More 
general moves can be derived involving both exchange 
and some modification of the exchanged states 



{Si,s 2 }->{s;(S2),s 2 (si)} 



(5) 



Such moves are allowed as long as they preserve detailed 
balance and - in combination with the independent sim- 
ulations - do not violate global ergodicity. The exponent 
of the acceptance probability for such generalized replica 
exchange moves is no longer given by the simple form in 
Eq. (0} but by 
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(6) 



Molecular dynamics, i.e. determining the time evolu- 
tion of a classical many-body system by numerically solv- 
ing the equations of motion, is an intrinsically determinis- 
tic approach. At first glance, this property seems to pre- 
clude its incorporation into the replica exchange scheme 
sketched above, since no stochastic sampling appears to 
be involved. However, already early on MD simulations 
were viewed also as a stochastic sampling procedure by 
considering the velocity field as a heat bath [T3, [H| : A 
continuous appropriate rescaling of the velocities leads to 
a correct canonical sampling of all properties depending 
only on the coordinates of the particles [l2|, EH] • In this 
way the system state is defined by the coordinates only, 
S = x, and the velocity field in a way acts as stochastic 
sampler. The above replica exchange scheme can then be 
applied straightforwardly. Care has to be taken, however, 
to properly treat the velocity field as stochastic sampler 
also in this situation, e.g. like in Ref. by randomly 
reinitializing the velocities upon exchange according to 
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(7) 



with N being the number of degrees of freedom. 

By considering the velocity field as a heat bath one 
gives up the idea of a trajectory. For this and other pjj] 
reasons, it is more common in canonical MD simulations 
to control the system by an additional thermostat algo- 
rithm [ToL ITU fT^ I . In such an implementation of canoni- 
cal MD the system state is given by coordinates and ve- 
locities together, S — (x, v), the thermostat acting as the 
external stochastic sampler. Consequently, in replica ex- 
change moves the full state of the system, i.e S = (x, v), 
has to be considered. The contributions to the energy 
can be separated as 



E(S) = E pot {x) + E kin (v) 



(8) 



As a consequence, the acceptance probability of Eq. ^ 
takes now on the form 



with 



ln(W) = ApAE pot + AK (9) 



AK = (fc - ft) [E kln (v 2 ) - E km ( Vl )} . (10) 



Since the kinetic energy reflects the simulation tempera- 
ture, see Eq. ([7]), in general AK will be negative. Hence, 
a naive application of Eq. (frj is hampered by a large 
detrimental contribution from the kinetic energy differ- 
ence yielding a very low acceptance probability, Eq. ([J). 



Moreover, accepted moves also lead to velocities that are 
not characteristic for the new temperatures, pushing the 
system out of equilibrium. 

One way to avoid such problems is to turn to general- 
ized exchange moves that control the possibly large fluc- 
tuations in the velocity fields. A re-scaling of the velocity 
fields in the move 

{(xi,vi),(x 2) v 2 )} -> {(x 2) r _1 V2),(xi,rvi)} , (11) 

leads now to a contribution of the kinetic energy to the 
acceptance probability of Eq. ([9]) that is given by 

AK(r) = (ft - ftr- 2 )(£ fcm (v 2 ) - r 2 £ fcm (vi)) . (12) 

where we have used the property E k i n (rv) = r 2 E kin {w). 
In order to optimize the acceptance of replica exchange 
moves one can adjust the scale r in such a way that the 
kinetic energy contribution to the Metropolis term van- 
ishes. The condition AK{r) — has two solutions. Fol- 
lowing Ref. [i| one can choose 



(13) 



This choice of r leads to ft — ftr - - = in Eq. (fT2f . but 
obviously this is not the only possibility. An alternative 
is to adjust r in a way that E k i n (y 2 ) — r 2 E k i n (y{) = 0. 
This can be realized by setting 




\/E U n{vi)lE kin {v x ) = y/fa/fi . (14) 



Here, T — 2kBE k i n is the instantaneous temperature of 
the velocity field, as opposed to the thermostat temper- 
ature T. Obviously, in the thermodynamic limit T — > T, 
and both scalings become identical. However, in finite 
systems T fluctuates around T and, in general, one has 
r E ^ rp. 

As both scalings preserve detailed balance [14|, they 
lead to equally valid but different replica exchange moves. 
Both scalings are optimal in the sense that they render 
acceptance rates independent of fluctuations in the veloc- 
ity field. Current implementations usually employ only 
a single variant. Choosing randomly among both moves 
can increase mixing in replica exchange, albeit without 
further increasing acceptance rates. 

Since AK(r) — > —00 for r — > as well as for r — > 00, 
and AK(rp) = = AK{te), there exists a regime of r 
where AK{r) > 0. It is therefore tempting to maximize 
AK(r) as this will increase the acceptance probability of 
an exchange move. From the condition 
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FIG. 2: Sketch of the coordinate part of the phase space 
for constant energy replica exchange simulations. Exchange 
moves are possible in the overlap region. 



which leads to a positive maximal contribution of the 
kinetic energy given by 



( E rC> /A B 2 

: t : 




FIG. 3: Dependence of the average acceptance probability 
{pai) on the distance of the energy shells, see Eq. (|23[1 . in 
microcanonical replica exchange simulations. 



(r op t) = (y/3i£fcm(vi) - V/3 2 ^n(v 2 )) 2 . (17) 



AK 



Hence, with such a re-scaling of velocities the fluctuations 
of the instantaneous temperatures of the velocity fields 
can be utilized to increase the acceptance probability of 
exchange moves, allowing for larger jumps in potential 
energy. Note, however, that the re-scaling according to 
Eq. (|16p preserves detailed balance only in the thermody- 
namic limit, i.e. where T — ^ T holds and AK(r opt ) = 
anyhow. For finite systems and during equilibration, i.e. 
where T deviates from T and AK(r opt ) > 0, detailed bal- 
ance is violated, albeit to a lesser degree the larger the 
system is. For this reason, the scaling of Eq. (|16p has to 
be used with care. 

The above theoretical framework allows us to intro- 
duce also an interesting new variant of replica exchange 
MD. Here the stochastic sampler is just the regular MD 
simulation without a heat bath. Provided the dynam- 
ics is ergodic, a microcanonical MD simulation can be 
viewed as a constant sampling on the energy surface, 
E = Ekin + Ep t = const. The corresponding weight 
function is w(S) = l/\fl(E)\, where Sl(E) denotes the 
phase space of the hypersurface of constant energy at 
E. It is well-known that microcanonical MD exhibits 
slow equilibration, and independently sampled trajecto- 
ries should be combined to ensure better statistics. How- 
ever, this approach can be extended readily into a replica 
exchange scheme for speeding up equilibration on several 
energy surfaces together. We assume E\ < E% in the 
following. The move set is a generalization of Eq. (Tl"Tj) 



involving two different rescaling factors r\ and r 2 



{(xi,vi), (x 2 ,v 2 )} -> {(x 2 ,r 2 v 2 ), (xi,nvi)} 



(18) 
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' E 2 ,i — Spot(xi. 2 ) 

£-1,2 — E po t{x\fl) 



(19) 



Such moves are possible for £ , pot (x 2 ) < Ex, see Fig. 2; 
Epotfax) < E2 automatically holds. We note that this re- 
striction does not violate detailed balance. Furthermore, 
the combination with the regular MD simulation ensures 
ergodicity. 

A fascinating aspect of this scheme is that the accep- 
tance probability for an allowed move is always one, since 
both weight functions are constant. No other scenario 
allows for such a high acceptance rate. Using reweight- 
ing techniques [17| , canonical properties can be obtained 
from simulations at several different energy values. More- 
over, constant energy surface simulations may be of in- 
terest in their own right [lOLIlll |. e.g. for comparison with 
recent molecular beam experiments 18 1. 

Practical acceptance probabilities will be somewhat 
smaller than one, (jpm) < 1, due to the region of for- 
bidden moves. Using the distribution of the potential 
energy on energy shell Ei, PEi{E pot ), the average accep- 
tance probability is given by 



(P 



dE pot PE 2 (Ep t) 



(20) 



For classical trajectories of total energy E equipartition 
of average kinetic and potential energy usually holds, 



(E k 
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(21) 
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with E, 



(mm) 
pot 



the energy of the lowest energy configura- 
tion. Assuming a Gaussian distribution for the potential 
energy, 



PEiiEpot) oc exp - (E pot - (E pot ) E J /A|. 



(22) 



we obtain the dependence of the average acceptance 
probability (pm) on the energy difference E 2 — E\ 



, . 1 , ( E 2 - El E,-El7 nr 
M-2 erfc ^ 



(23) 



shown in Fig. 3. In particular, (pm) will be at least one 
half if the average potential energy at E 2 is smaller than 
Ei, {Epot) E ^ < E\. This criterion is equivalent to the 
energy difference being at most equal to twice the average 
kinetic energy at E 1 , E 2 -E x < Ei-E [ p ™ n) = 2 (E kin ) Ei . 
Figure 3 shows that (pm) will rapidly approach one upon 
decreasing energy difference. 

This last example demonstrates the wide applicabil- 
ity of generalized replica exchange move sets. It also 
demonstrates the striking advantage of replica exchange 
over earlier approaches like simulated tempering [19| . In 
the latter, additional parameters reflecting free energy 
differences are of utmost importance. Their determina- 
tion is tedious and approximations [2fJ are useful only 



in certain limiting cases. In replica exchange such nor- 
malization constants simply drop out due to the form of 
the acceptance probability, Eq. ([3]). Constant energy sur- 
face simulations as sketched above could be approached 
in simulated tempering only with a solid knowledge of 
the phase space ratios, |J7(£ , 2)|/|51(i?i)|. 

In summary, we have sketched the statistical physics 
framework of applying the replica exchange technique to 
MD simulations. Generalized move sets, in particular 
appropriate rescaling of the velocity fields, allow opti- 
mization of the acceptance probability as well as new ap- 
proaches. Together with an optimization of the temper- 
ature spacing to increase replica flow (l5l . Il6j optimized 
acceptance probabilities will lead to shorter simulation 
times in canonical replica exchange MD simulations. Mi- 
crocanonical replica exchange MD simulations are intrin- 
sically optimized and will provide new insights. 
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